# Alex Gazmararian
# agazmararian@gmail.com
source("analysis/visibility/analysis_config.R")

g <- readRDS(here("data", "output", "visibility_analysis.rds"))


# Use treatment variables from config, but remove _q suffix to get continuous variables
treat_var_re <- gsub("_q$", "", treat.vars[1])  # Renewable energy treatment variable
treat_var_mfg <- gsub("_q$", "", treat.vars[2]) # Manufacturing treatment variable

resid.red <- feols(reformulate(covar.list, response = treat_var_re), data = g, fixef = c("state")) %>%
  resid() %>%
  data.frame(value = .)

resid.mfg <- feols(reformulate(covar.list, response = treat_var_mfg), data = g, fixef = c("state")) %>%
  resid() %>%
  data.frame(value = .)

# Which states have the most variation
state.df <- data.frame("state" = g$state, "re" = resid.red$value, "mfg" = resid.mfg$value)
state.df %>%
  filter(re > 10 | re < -10) %>%
  group_by(state) %>%
  tally() %>%
  arrange(desc(n))
state.df %>%
  filter(mfg > 10 | mfg < -10) %>%
  group_by(state) %>%
  tally() %>%
  arrange(desc(n))

x_min <- min(resid.red$value, resid.mfg$value, na.rm = TRUE)
x_max <- max(resid.red$value, resid.mfg$value, na.rm = TRUE)
y_max <- 2000

p.dist.re <- resid.red %>%
  filter(!is.na(value)) %>%
  ggplot(aes(x = value)) +
  geom_histogram(color = "black", bins = nclass.scott(resid.red$value)) +
  theme_classic(base_size = 10) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, y_max)) +
  scale_x_continuous(expand = c(0, 0), limits = c(x_min, x_max)) +
  annotate("text", x = 50, y = 1000, label = paste0("Within-state SD=", round(sd(resid.red$value),1)), hjust = 0) +
  labs(x = "Within-state residualized distance to project (km)", y = "Count",
  title = "Renewable energy")

p.dist.mfg <- resid.mfg %>%
  filter(!is.na(value)) %>%
  ggplot(aes(x = value)) +
  geom_histogram(color = "black", bins = nclass.scott(resid.mfg$value)) +
  theme_classic(base_size = 10) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, y_max)) +
  scale_x_continuous(expand = c(0, 0), limits = c(x_min, x_max)) +
  annotate("text", x = 50, y = 1000, label = paste0("Within-state SD=", round(sd(resid.mfg$value),1)), hjust = 0) +
  labs(x = "Within-state residualized distance to project (km)", y = "Count",
  title = "Green manufacturing")

p.dist.out <- p.dist.re + p.dist.mfg + 
  plot_layout(ncol = 2) +
  plot_annotation(tag_levels = "A")

save_pnas_pdf(
    p.dist.out,
    here("output", "pnas", "figures", "fig_S1_dist_respondent.pdf"),
    width = "double",
    height = 6
    )
message("Saved figure to ", here("output", "pnas", "figures", "fig_S1_dist_respondent.pdf"))
